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Capillary filling using Lattice Boltzmann Equations: the case of 
multi-phase flows 



o 
o 

(N 



o 

u 

d 



> 
in 

On 
P 

o 
o 



X 



F. Diotallevi^, L. Biferale^, S. Chibbaro^, A. Lamura^, G. Pontrelli^, M. Sbragaglia^, S. Succi^, and F. Toschi^ 

^ Istituto per le Applicazioni del Calcolo CNR, Viale del Policlinico 137, 00161 Roma, Italy. 

^ Dept. of Physics and INFN, University of Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma, Italy. 

^ Dept. of Mechanical Engineering, University of Tor Vergata, Viale Politecnico 8, Rome, Italy. 

^ Istituto per le Applicazioni del Calcolo CNR, Via Amendola 122/d 70126 Bari, Italy. 

^ Department of Applied Physics, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. 
Received: date / Revised version: date 

Abstract. We present a systematic study of capillary filling for multi-phase flows by using mesoscopic 
lattice Boltzmann models describing a diffusive interface moving at a given contact angle with respect to 
the walls. We compare the numerical results at changing the density ratio between liquid and gas phases, 
Sp/p and the ratio, S^/H, between the typical size of the capillary, if, and the interface width, S^. It is 
shown that numerical results yield quantitative agreement with the Washburn law when both ratios are 
large, i./e. as the hydrodynamic limit of a infinitely thin interface is approached. We also show that in 
the initial stage of the filling process, transient behaviour induced by inertial effects and "vena contracta" 
mechanisms, may induce significant departure from the Washburn law. Both effects are under control in 
our lattice Boltzmann equation and in good agreement with the phenomenology of capillary filling. 
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1 Introduction 

The physics of capillary filling is an old problem, originat- 
ing with the pioneering works of Washburn [l] and Lucas 
[2]. Recently, with the explosion of theoretical, experimen- 
tal and numerical works on microphysics and nanophysics, 
the problem attracted a renewed interest [3l[4l[5l[6] . Capil- 
lary filling is a typical "contact line" problem, where the 
subtle non-hydrodynamic effects taking place at the con- 
tact point between liquid-gas and solid phase allows the 
interface to move, pulled by capillary forces and contrasted 
by viscous forces. Usually, only the late asymptotic stage 
is studied, leading to the well-known Washburn law, which 
predicts the following relation for the position of the in- 
terface inside the capillary: 



capillary time tcap 
law: 



Hfi/j. This leads to the universal 



cos{0) 



(2) 



-fHcos{0)~ 
3/x 



(1) 



where 7 is the surface tension between liquid and gas, 
is the static contact angle, fi is the liquid viscosity, H 
is the channel height and the factor 3 depends on the 
geometry of the channel (here a two dimensional geometry 
given by two infinite parallel plates separated by a distance 
H - see fig. [1]). The above expression can be recasted in 
dimensionless variables t = t/tcap and z = z/H^ being the 
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As already remarked in many works [5], the asymptotic 
behaviour ([2]) is obtained under the assumption that (i) 
the inertial terms in the Navier- Stokes evolution are neg- 
ligible, (ii) the instantaneous bulk profile is given by the 
Poiseuille flow, (iii) the microscopic slip mechanism which 
allow for the movement of the interface is not relevant 
for bulk quantities (as the overall position of the inter- 
face inside the channel), (iv) inlet and outlet phenomena 
can be neglected (limit of infinitely long channels); (v) the 
liquid is filling in a capillary, either empty or filled with 
gas whose total mass is negligible with respect to the liq- 
uid one. In the following, we will address all these effects 
and show to which extent they can described by using a 
mesoscopic model for multiphase flows based on the dis- 
cretized version of Boltzmann Equations in a lattice. The 
model here used is a suitable adaptation of the Shan-Chen 
pseudo-potential LBE [7J with hydrophobic/hydrophilic 
boundaries conditions, as developed in Other mod- 

els with different boundary conditions and/or non-ideal 
interactions have been also used in [TO]. 
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2 LBE for capillary filling 

The geometry we are going to investigate is depicted in 
fig. ([T|). The bottom and top surface is coated only in 
the rigth half of the channel with a boundary condition 
imposing a given static contact angle [8j ; in the left half we 
impose periodic boundary conditions at top and bottom 
surface in order to have a flat liquid-gas interface which 
should mimic an "infinite reservoir". Periodic boundary 
conditions are also imposed at the two lateral sides such 
as to ensure total conservation of mass inside the system. 



2.1 LBE algorithm for multi-phase flows 

We start from the usual lattice Boltzmann equation with 
a single-time relaxation pT| [T2]: 
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At 

TB 



{fi{x,t) 



(3) 

where fi(x^t) is the kinetic probability density function 
associated with a mesoscopic velocity c^, tb is a mean 
collision time (with At a time lapse), //^^^ (p, pu) the equi- 
librium distribution, corresponding to the Maxwellian dis- 
tribution in the continuum limit. From the kinetic distri- 
butions we can define macroscopic density and momentum 
fields as [TTIO: 



P{^) =^ fi{^)'^ pu{x) = ^cifi{x). (4) 



For technical details and numerical simulations we shall 
refer to the nine-speed, two-dimensional 2DQ9 model [llj. 
The equilibrium distribution in the lattice Boltzmann equa- 
tions is obtained via a low Mach number expansion of the 
equilibrium Maxwellian [ll,12j. In order to study non- 
ideal effects we need to supplement the previous descrip- 
tion with an interparticle forcing. This is done by adding 
a suitable Fi in ([3]). In the original model [7J, the bulk in- 
terparticle interaction is proportional to a free parameter 
(the ratio of potential to thermal energy), ^5, entering the 
equation for the momentum balance: 

Fi = -gbclY,^i\''i\^M^^t)^i^ ^ ciAt,t)ci (5) 



being w{\ci\'^) the static weights for the standard case of 
2DQ9 [llj and i^{x^t) = ip{p{x^t)) the pseudo-potential 
function which describes the fiuid-fluid interactions trig- 
gered by inhomogeneities of the density profile (see [TllHl 
[9] for details). 

One may show [8^9j that the above pseudo-potential, 
leads to a non-ideal pressure tensor given by (upon Taylor 
expanding the forcing term): 
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Fig. 1. Geometrical set-up of the numerical LBE. The 2 di- 
mensional geometry, with length 2L and width H, is divided in 
two parts. The left part has top and bottom periodic bound- 
ary conditions such as to support a perfectly flat gas-liquid 
interface, mimicking a "infinite reservoir" . In the right half, of 
length L, there is the true capillary: the top and bottom bound- 
ary conditions are those of a solid wall, with a given contact 
angle [8 . Periodic boundary conditions are also imposed at 
the west and east sides. 



where Cg is the sound speed. This approach allows the def- 
inition of a static contact angle 6>, by introducing at the 
walls a suitable value for the pseudo-potential ip{pw) [H], 
which can span the range e [0^ : 180^]. Moreover, it 
also defines a specific value for the surface tension, 7/^, 
via the usual integration of the offset between normal and 
transverse components of the pressure tensor along the 
liquid-gas interface allows for[7^8^9j. 

As to the boundary conditions on the Boltzmann popu- 
lations, the standard bounce-back rule is imposed. One 
can show that the bounce-back rule gives no-slip bound- 
ary conditions up to second order in the Knudsen num- 
ber in the hydrodynamical limit of single phase flows [13j. 
In presence of strong density variation, close to the walls 
across the interface, the velocity parallel to the wall may 
develop a small slip length (of the order of the interface 
thickness, Ag ex J^) which in turn, allows for the interface 
to move. It is difficult to control exactly the phenomenon, 
because even imposing an exact no-slip boundary condi- 
tions at the wall [14J, the model will develop non triv- 
ial dynamics at the first node away from the wall, where 
both condensation/evaporation phenomena and/or spuri- 
ous currents may conspire, leading to an overall non-zero 
slip velocity. For the scope of controlling the capillary fill- 
ing, one may reabsorb all these effects within the usual 
Maxwell slip boundary conditions: Us = XsdnU- It is easy 
to show that in presence of a slip velocity, the Poiseuille 
profile becomes: 



u{y) = 6 



y{H -y)^XH 



6Xs/H 



(7) 



where the velocity of the front must be identified with 
the mean velocity, u = 1/H dyu{y) = z. Therefore, 
Washburn law ^ becomes: 



(6) 



^(1 + 6^). 
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2.2 Corrections to the Wahsburn law 



As already remarked many yeas ago [15j, the Wahsburn 
law ([2j) can be valid only if inertial forces are negligible 
with respect to the viscous and capillary ones. This cannot 
be true at the beginning of the filling process, where strong 
acceleration drives the interface inside the capillary. How- 
ever, putting reasonable numbers for a typical microde- 
vices experiments with water {H 1/im, 7 0.072A/'/m, 
pi 10~^kg/m^^ ji l{)~^Ns/m?)^ one realizes that the 
transient time, Tdiff = H^pi/ ji^^ is usually very small, of 
the order of a few nanoseconds, and therefore negligible 
for most experimental purposes. Another important effect 
which must be kept in mind when trying to simulate cap- 
illary filling, is the unavoidable "resistance" of the gas oc- 
cupying the capillary during the liquid invasion. This is a 
particular "sensitive" question, because reaching the typi- 
cal 1 : 1000 density ratio between liquid and gas of exper- 
imental set up represents a challenge for most numerical 
methods, particularly for multiphase Lattice Boltzmann 
which typically operates with density ratios of the order 
1 : 10 or 1 : 100. In order to take in to account both ef- 
fects, inertia and gas dynamics, one may write down the 
balance between the total momentum change inside the 
capillary and the force acting on the liquid+gas system: 



d{zM{t)) 
Jt 



(9) 



where M{t) = Mg + Mi is the total mass of liquid and gas 
inside the capillary at any given time. The two forces in the 
right hand side correspond to the capillary force, Fcap = 



2^cos{6)^ and to the viscous drag F^. 



-2{pg{L-z)^ 



jiiiz)dnu{0). Following the notation of fig. ([I]) and the ex- 
pression for the velocity profile ([7j) one obtains the final 
expression (see also [16] for a similar derivation, without 
considering the slip velocity): 



{pg{L - z) ^ piz)z ^ {pi 
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In the above equation for the front dynamics, the terms 
in the LHS take into account the fluid inertia. Being pro- 
portional either to the acceleration or to the squared ve- 
locity, they become negligible for long times. Washburn 
law plus the slip correction (|8|) is therefore correctly re- 
covered asymptotically, for t ^ 00, and in the limit when 
pg/pi 0. The above equation is exact, in the case where 
evaporation-condensation effects are negligible, i.e. when 
the gas is pushed out of the capillary without any interac- 
tion with the liquid. This is not always the case for most 
of the mesoscopic models available in the literature [7^14J, 
based on a diffusive interface dynamics [18j. As we will see, 
only when either the limit of thin interface S^/H is 
reached or when the gas phase is negligible, Sp/p — > 00, the 
dynamics given by (p!Q|) is correctly recovered. Otherwise, 
deviations are observed, which are induced by condensa- 
tion/evaporation effects, which may result in significant 
departure from the Poiseuille profile inside the gas phase. 




Fig. 2. Evolution of the front in adimensional variables. Po- 
sition of the front, z{t), versus t for H / {Ax) — 15,61,121 
and L/{Ax) = 1200,1500,2000 respectively, and with 5(0) = 
10Z\x. The solid curve is the theoretical solution obtained by 
integration of ea. ()10|) with \s — 2. The data have been ob- 
tained with = 5 which corresponds to pg = 0.157 and 
pi — 1.92 (in LBE units). Let us notice that the solutions 
of eallOldoes not show any sensitive variations on H for those 
times and distances here explored. In the inset we show the 
position of the front for (H/Ax) — 15,61, 121 (same symbols) 
normalized with the adimensional asymptotic solution of eallOl 

given by expression 0: Zasym{t) (^i ^ 6^)t, notice 

that the departure from the predicted asymptotic Washburn 
law is never larger then %10 even for small channel width, and 
becomes almost negligible already for H — 121 Ax 



2.3 Numerical Results 

In fig. ([2|) we show the behaviour of the front position, 
z(t), as a function of time for a given contact angle {0 = 
55^ ±3^), a given density ration Sp/p =11 and a given 
surface tension 7 = 0.0569 (in LBE units), at varying the 
channel width from H = 15Ax up to H = 121Ax. 
As one can see the numerical results tends to be in good 
agreement with the solutions of (p!Q|) , only for large enough 
values of i.e. only when the interface becomes thin 
enough. For small to moderate values of the overall 
asymptotic trend is only qualitatively reproduced by the 
Washburn law (plus slip effects) (|8]) with deviations which 
may be of the order of 10% in the prefactor (see inset of 
the same figure). 

Similarly, increasing the surface tension and the dp/ p 
factor, leads to a early convergence towards the asymp- 
totic Washburn law and to the solution of (p!Q|) even for 
small channel width H. This is shown in fig. ([3]) where, at 
fixed H = 31Z\x, we increase 5p/ p and the Washburn law 
is approached better and better. 

From both figure ([2][3]) one may notice that for short 
filling time, t < Tdiff ^ strong deviations from the Wash- 
burn law are detected, even considering the extra effects 
induced by the inertial terms of (p!Q|) . This slowing down at 
the early stage of the filling is indeed mainly induced by a 
sort of vena contracta term [17 , reflecting the non-trivial 
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-- Eq. (10) with 6p/p=ll 






Eq. (10) with 8p/p=34 
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Eq. (10) with "Vena Contracta" 
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Fig. 3. Same plot of fig. {[2]) but for fixed H = 31Ax at chang- 
ing the density ration. We have for Qb = 5, Sp/p — 11 and 
for Qh — 6, 8p/ p — 34. Notice that the second case is already 
enough to have a very good agreement with the solution of {[TO)) 
(adimensionalized as explained in the text) also at this small 
channel width (solid line = 6, dashed line = 5). 



matching between the reservoir and the capillary dynam- 
ics at the inlet. This term, has been argued to be describ- 
able by an additive apparent-mass correction, cpiHz, to 
the LHS of (dni). 

In Fig. (|4]) we show an enlargement of fig ([2]) for small 
filling time superposed with the results of a numerical in- 
tegration of eq. (p!Q|) with a phenomenological value c = 30 
for the vena contracta factor. As one can see, the agree- 
ment between the numerics and the evolution of ([TQ|) is 
now excellent also at short times. 



3 Conclusions 



The present study shows that Lattice Boltzmann mod- 
els with pseudo-potential energy interactions are capa- 
ble of reproducing the basic features of capillary filling, 
as described within the Washburn approximation. Two 
conditions for quantitative agreement have been identi- 
fied: i) a sufficiently high density contrast between the 
dense/light phase, pi/pg > 10 and a sufficiently thin in- 
terface, 5£^/H < 0.1. Both conditions can be met within 
the current LB methodology, although it would clearly 
be desirable to extend the LB scheme in such a way to 
achieve density contrasts in the order of 1 : 1000 (the cur- 
rent state-of-the-art is approximately 1 : 50) and interface 
widths of the order of the lattice spacing Ax (current val- 
ues are about 5 — 10 Ax). 

The present results set the stage for future computa- 
tional studies aimed at identifying optimal interface func- 
tionalization strategies, based on physical, chemical and 
geometrical coating processes. Work along these lines is 
currently underway. 



Fig. 4. Enlargement of the early stage evolution during the 
filling process for H = 121 Ax (in LBE units), with the vena 
contracta term choosing an optimal value for c = 30. The front 
position and time are make dimensionless normalizing with tcap 
and with H. 
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